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ABSTRACT 

A  methodology  for  quantitative,  directional  validation  of  a  long-term  wave  model  hindcasl  is  described 
and  applied.  Buoy  observations  are  used  as  ground  iruifi  and  the  method  does  not  require  the  application 
of  a  parametric  model  or  data-adaptive  method  to  ihe  observations.  Four  frequency  ranges,  relative  to  the 
peak  frequency,  are  considered.  The  validation  of  the  hindcasl  does  not  suggest  any  systematic  bias  in 
predictions  of  directional  spreading  at  or  above  the  spectral  peak.  Idealized  simulations  are  presented  to  aid 
in  the  interpretation  of  results. 


1.  Introduction 

a.  Background 

1 )  iMPORTANCE/RELE VANCE 

Principal  wave  direction,  quantified  as  a  mean  nr 
peak  value,  is  of  obvious  importance  to  wave  predic¬ 
tion.  Directional  distribution  about  the  mean  or  peak 
direction  is  also  very  important  for  wave  modeling,  It 
can  have  a  large  impact  on  the  prediction  of  swells, 
since  il  determines  how  far  and  wide  the  swells  will 
disperse.  Nonlinear  interactions  computed  by  a  wave 
model  are  sensitive  to  the  directional  distribution  of 
energy.  Further,  as  wave  model  dissipation  terms  with 
more  sophisticated  directional  dependency  are  devel¬ 
oped,  we  can  expect  that  directional  spreading  will  have 
greater  influence  on  the  modeled  source  term  balance 
and,  thus,  total  energy. 

2)  Present  capability 

Validations  of  modeled  peak  or  mean  wave  direction 
in  the  Literature  typically  show  good  skill,  though  the 
response  of  a  third-generation  wave  model  to  rapidly 
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turning  winds  is  a  concern  (e,g..  Young  el  al.  1987).  The 
ability  of  third-generation  models  to  accurately  predict 
the  width  of  the  directional  distribution  is  poorly  un¬ 
derstood.  Indeed,  as  is  described  in  a  companion  manu¬ 
script  (Rogers  and  Wang  200h,  hereafter  RW),  evalua¬ 
tions  in  the  literature  show  very  little  consensus. 

h .  Model  description 

The  so-called  third-generation  (3G)  of  spectral  wave 
models  calculate  wave  spectra  without  a  priori  assump¬ 
tions  regarding  spectral  shape.  For  this  investigation, 
we  use  the  Simulating  Waves  Nearshore  model 
(SWAN:  Booij  et  al.  1999).  SWAN  is  a  3G  model  de¬ 
signed  to  address  the  excessive  computational  expense 
of  applying  predecessor  3G  models  (such  as  WAM: 
WAMDI  Group  1988}  in  coastal  regions.  The  govern¬ 
ing  equation  of  SWAN  and  most  other  3G  wave  models 
is  the  action  balance  equation.  In  Cartesian  coordi¬ 
nates,  the  action  balance  equation  is 

BN  flC^/V  BC  VN  aC-JV  S 

— -  +  ^ — —  +  — - —  +  — r~ —  +  — — —  = 

dt  dx  dy  Ba  36  it 

where  a  is  the  relative  frequency,  which  is  the  wave 
frequency  measured  from  a  frame  of  reference  moving 
with  a  current,  if  a  current  exists:  N  is  the  wave  action 
density,  equal  to  the  energy  density  divided  by  the  rela¬ 
tive  frequency  (N  —  E/tr):  6  is  the  wave  direction;  CK  is 
the  wave  action  propagation  speed  in  (a;  y,  ft)  space: 
and  S  is  ihe  total  of  source-sink  terms  expressed  as 
wave  energy  density.  In  deep  water,  lire  right-hand  side 
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of  the  governing  equation  is  dominated  by  three  terms: 
S  Sin  +  Sn |  4-  .\|K  (input  by  wind,  four  wave  nonlinear 
interactions,  and  dissipation,  respectively).  These  three 
deepwater  source-sink  terms  are  discussed  at  several 
points  later  in  this  manuscript,  SWAN  also  includes 
physical  processes  associated  with  intermediate-depth 
and  shallow  water  (e.g,.  bottom  friction,  depth- limited 
breaking). 

c.  Objective 

It  has  become  increasingly  common  fora  wave  mod¬ 
eler  to  have  at  his  or  her  disposal  directional  wave  ob¬ 
servations  within  a  model  computational  domain.  This 
often  leads  to  an  expectation — perhaps  a  naive  expec¬ 
tation — that  the  wave  modeler  can  readily  use  these 
observations  to  validate  the  model.  Unfortunately,  vali¬ 
dating  a  model  using  directional  observations  is  much 
less  straightforward  than  traditional  validations  of  wave 
height  or  peak  period.  For  example,  what  if  the  model 
in  question  is  a  long-term  simulation  with  continuous 
directional  observations?  How  could  one  perform  a 
meaningful  validation  that  is  compact  enough  to  be  pre¬ 
sented  to  others?  How  far  can  one  go  in  condensing 
these  comparisons?  At  what  point  do  the  comparisons 
become  meaningless  or  misleading? 

This  study  was  initiated  with  three  general  objectives: 
I)  to  review  the  history  and  state  of  the  art  for  direc¬ 
tional  wave  validation  methods,  2)  to  design  a  valida¬ 
tion  methodology/strategy  best  suited  for  a  specific 
model  application,  and  3)  to  characterize  model  behav¬ 
ior  in  that  specific  application.  The  first  general  objec¬ 
tive  is  addressed  only  briefly  in  this  manuscript:  a  com¬ 
panion  manuscript  (RW)  provides  a  more  detailed  re¬ 
view.  The  specifics  of  the  second  and  third  general 
objectives  are  given  here. 

1)  Requirements  on  vai  idation  design 

The  major  challenge  of  this  study  is  in  the  design  of 
a  validation  method.  Since  we  are  allowed  the  luxury 
here  of  focusing  primarily  on  directional  wave  valida¬ 
tion,  we  were  not  satisfied  with  the  simplest  and  most 
obvious  method,  which  is  to  use  buoy  data  and  a  para¬ 
metric  model  or  data-adaptive  method  to  create  direc¬ 
tional  spectra,  and  make  qualitative  side-by-side  com¬ 
parisons  with  model  directional  spectra  at  a  few  instants 
in  time.  Rather,  we  have  fairly  specific  self-imposed 
requirements  on  the  validation  method. 

The  first  requirement  is  that  it  be  a  long-term  direc¬ 
tional  validation,  without  extensive  manipulation  of  the 
output,  for  example  to  isolate  the  pure  vvindsca  events. 
Usually,  when  directional  spreading  is  a  primary  focus, 
the  investigators  focus  on  specific  events.  This  leads  to 
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uncertainties  with  regard  to  generality:  is  a  conclusion 
specific  to  the  event,  or  is  it  a  systematic  symptom  of 
the  model  physics?  We  address  this  limitation  using  a 
relatively  long  simulation. 

The  second  requirement  is  to  develop  and  employ  a 
method  of  evaluation  of  model  directional  skill  that  is 
t tuaniitative*  in  other  words  a  comparison  of  model  and 
observation  value  pairings  from  which  statistics  may  be 
calculated:  traditionally,  this  is  presented  as  a  scatter- 
plot  comparison  of  modeled  and  observed  values.  Since 
many  comparisons  in  the  literature  are  short-term  com¬ 
parisons.  it  is  possible  to  simply  present  modeled  and 
observed  two-dimensional  spectra  side  by  side,  thereby 
avoiding  the  necessity  of  condensing  results.  With  long¬ 
term  simulations,  it  is  necessary  to  condense  results 
somehow. 

The  third  requirement  is  that  the  method  of  evalua¬ 
tion  of  model  directional  skill  also  utilizes  observational 
data  as  they  are  given t  rather  than  applying  a  parametric 
model  |e.g*.  the  well-known  con2''  model:  Longuet- 
Higgins  el  al.  1963);  Cartwright  (1963)]  or  a  data- 
adaptive  method.  Two  popular  methods  are  the  maxi¬ 
mum  likelihood  estimator  [referred  to  as  VILE  or 
MLM;  Capon  cl  al.  (1%7)  and  Oltman-Shay  and  Guza 
(1984)]  and  the  maximum  entropy  method  (MEM; 
Lygre  and  Krogstad  1986)  to  transform  the  observa¬ 
tional  data  into  a  subjective  directional  spectrum.  For 
discussion  and  a  description  of  this  subjectivity,  we  re¬ 
fer  the  reader  to  Kuik  el  al.  (1988)  and  Benoit  el  al. 
(1997), 

Our  fourth  requirement  is  that  the  observational  data 
be  taken  from  a  buoy.  Other  data  sources,  such  as  ra¬ 
dar.  have  been  used  with  success  in  the  past  for  direc¬ 
tional  validation,  but  these  datasets  tend  to  have  more 
limited  availability  or  accessibility. 

Our  fifth  requirement  is  that  the  frequency  variation 
of  directional  spreading  be  considered,  as  opposed  to  a 
quantity  integrated  from  the  entire  frequency  range, 

2)  Model  characterization 

The  objective  is  to  quantitatively  determine  whether 
a  typical  3G  model  (SWAN;  Booij  et  al.  1999),  in  a 
typical  implementation,  has  a  systematic  tendency  to 
overpredict  or  underpredict  directional  spreading.  The 
Discrete  Interaction  Approximation  (D1A)  for  four- 
wave  nonlinear  interactions,  Sn |4,  is  the  approximation 
used  by  all  operational  3G  wave  models  today.  It  is  has 
been  demonstrated  a  number  of  times  in  the  literature 
that  this  approximation  leads  to  broader  directional 
spreading  than  would  be  obtained  with  more  rigorous 
calculations  (1  lasselmann  el  al,  1985;  Young  et  al.  1987; 
Young  and  Van  Vledder  1993;  Cardone  and  Resio 
1998;  Forristall  and  Ewans  1998,  etc,).  This  can  result  in 
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Table  L  Summary  of  notation. 


/ 


tr 

H 

<M/) 

N 

E 

F 

D{B) 

/,  and  /, 

U  /) 

%{f) 

nq(/) 

tf|.  /»,.  a?.  b> 


Frequency,  T  ■ 

The  relative  (intrinsic)  radial  frequency.  2t tT  1 
Direction  of  wave  propagation 

The  rms  circular  spreading  (note,  tr  here  is  unrelated  to  frequency) 

Mean  rrm  circular  spreading:  "mean”  here  refers  to  some  integration  over  frequencies 
Two-dimensional  spectral  wave  action  density,  N{  f.  it) 

Two-dimensional  spectral  variance  density.  Fi  t,  0) 

One -dimensional  spectral  variance  density.  F(  f ) 

Dimensionless  directional  distribution  at  a  particular  frequency:  integrates  to  unity 
Lower  and  upper  bounds  of  a  frequency  integration 

Mean  wave  direction,  taken  as  the  circular  centroid  of  Pp9,  and  denoted  o,( f)  in  ND1K  notation 

Peak  wave  direction,  the  peak  of  D{0):  generally  not  known,  except  in  the  context  of  a  model  of  some  sort 

Mean-mean  wave  direction,  which  has  been  integrated  across  some  frequency  range 

Parameter  related  to  directional  spreading,  denoted  r,  in  NDBC  notation 

Fourier  coefficients 


an  expectation  that  3G  wave  models  systematically 
overpredict  directional  spreading.  This  is  sometimes 
observed  in  the  literature,  hut  the  reverse  has  also  been 
found  (Jensen  et  a  I,  1995).  One  wonders  how  much  this 
"expectation"  has  influenced  prior  comparisons.  Long¬ 
term  comparisons  can  he  used  to  convincingly  argue  for 
or  against  this  broadening  effect.  We  know  from  die 
literature  that  in  cases  of  pure  windsea,  directional 
spreading  lends  to  follow  a  fairly  consistent  pattern 
relative  to  the  peak  frequency:  directional  spreading  at 
the  peak  is  relatively  narrow,  and  spreading  is  broader 
farther  from  (higher  and  lower  than)  the  peak.  A  sec¬ 
ondary  objective  is  to  verify  that  an  operationally  used 
wave  model  (SWAN  with  the  DIA  approximation  for 
four-wave  interactions)  adequately  reproduces  this  pat¬ 
tern  in  directional  spreading. 

3)  Prior  works 

A  number  of  methods  for  directional  validation  of 
wave  predictions  have  been  applied  over  the  years.  Due 
to  page  limits,  a  detailed  review  of  this  prior  work  is 
described  separately  (RW).  There  have  been  no  prior 
works  that  fit  the  five  requirements  described  above. 
There  have,  of  course,  been  a  number  of  studies  that 
share  some  similarities.  For  example,  Komen  et  al. 
(1994,  chapter  V.4).  Khandekar  et  al.  (1994),  and  For- 
ristall  and  Greenwood  (1998)  describe  the  validation  of 
directional  spreading  of  hindeasts  of  medium  ( 15  days) 
or  longer  duration.  Jensen  et  al.  ( 1995),  Forristall  and 
Greenwood  ( 1998),  Ardhutn  et  at.  (2003),  and  Wyatt  et 
al.  (2003)  all  include  quantitative  non-da  la-adaptive 
comparisons  lor  validation  of  hindcast  directional 
spreading  with  in  situ  data  as  ground  truth.  Forristall  el 
al.  (1978),  Komen  et  al.  (1984),  Toiman  (1991),  Forri- 
stall  and  Ewans  (1998),  Forristall  and  Greenwood 
( 1998),  and  Alves  and  Banner  (2003)  include  quantita¬ 


tive  validations  of  directional  spreading,  with  some  de¬ 
scription  of  the  variation  with  frequency.  Of  these,  only 
Forristall  et  a!.  (1978)  was  a  hindcast  of  a  specific  event, 
as  opposed  to  an  idealized  simulation,  and  showed  the 
frequency  variation  simply  by  choosing  a  few  specific 
instances  in  time.  For  further  descriptions  of  prior  work 
on  the  subjects  of  directional  metrics,  directional  model 
validation,  and  parametric  directional  distributions,  we 
refer  the  reader  to  RW. 

it.  Terminology 

The  two-dimensional  energy  density  spectrum  is 
defined  as  E(f.  9)  =  D(fm  0)F(/),  where  />(/,  9)  ts 
the  normalized  directional  distribution  and  F(  f )  is 
the  one- dimensional  energy  density  spectrum. 
The  function  />(/,  9)  is  normalized  such  that 

"Directional  spreading"  refers  to  the  degree  to  which 
a  directional  distribution  of  wave  energy  is  "broad.”  It 
does  not  refer  to  the  normalized  directional  distribution 
itself,  which  is  sometimes  referred  to  as  the  "directional 
spreading  function."  Notations  used  herein  are  given  in 
Table  1. 

e.  Organization  of  manuscript 

In  section  2,  the  methodology  of  this  study  (general 
validation  strategy  and  definition  of  metrics  used)  is 
described.  In  section  3,  an  idealized  case  is  examined  to 
isolate  the  effect  of  the  inaccuracy  of  the  Discrete  Inter- 
Action  approximation  for  four-wave  nonlinear  interac¬ 
tions.  In  section  4,  an  example  directional  validation  is 
presented  fora  hindcast  with  the  SWAN  model  in  Lake 
Michigan  during  fall  2002.  Results  are  summarized  in 
section  5,  Discussion  is  given  in  section  6  and  conclu¬ 
sions  in  section  7. 
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2.  Method 

a.  General  'strategy 

1 )  Ground  truth 

Buoy  data  are  the  “ground  truth”  of  this  study — 
specifically*  National  Data  Buoy  Center  (NDBC)  buoy 
45007  in  Lake  Michigan.  Directional  buoys  are  often 
the  most  cost-effective  method  of  obtaining  directional 
data  outside  the  surf  zone.  [In  depths  shallower  than 
around  150  m.  three-element  pressure  gauge  arrays  and 
p-u-v  gauges  can  be  cost-effective  methods  of  obtain¬ 
ing  information  that  is  essentially  the  same  as  that  from 
a  heave-pitch -roll  buoy.  Additional  elements  in  a  pres¬ 
sure  gauge  or  wave  staff  array  will  yield  higher- 
resolution  directional  data:  see  Young  (1994).] 

2)  N  ON  D  t  REtT  I  ON  A 1.  ACCURACY 

We  aim  to  perform  a  model  validation  in  which  di¬ 
rectional  characteristics  are  the  primary  focus.  Usually, 
when  directional  metrics  are  used  in  validation,  ihey 
are  secondary,  with  the  primary  focus  being  wave 
height,  wave  period,  and — more  rarely — frequency 
spectra.  Here,  we  want  to  devote  most  of  our  attention 
to  the  directional  issue.  We  do  this  by  taking  a  modeling 
system  that  has  consistently  good  skill  with  regard  to 
nundirectional  metrics. 

3)  Challenge:  Problem  complexity 

Our  objective  is  to  determine  the  feasibility  of  con¬ 
ducting  a  quantitative  directional  validation  of  a  long¬ 
term  hindcast.  Anticipating  that  is  a  major  challenge 
even  under  the  most  favorable  circumstances:  we  sim¬ 
plify  our  case  study  by  taking  the  following  steps. 

1 )  We  use  a  lake  as  our  lest  basin  (Lake  Michigan,  Fig. 
1);  thus,  the  wave  climate  is  dominated  by  windsea. 
Mixed  sea/swell  states  (identifiable  as  having  mul¬ 
tiple  peaks)  do  occur  (especially  when  the  wind 
shifts  rapidly),  hul  are  uncommon.  Certainly,  old 
swells  do  not  occur. 

2)  We  use  a  model  (SWAN)  that  has  proven  to  be 
skillful  in  predicting  nondirectional  spectra  at  this 
scale,  in  wind  sea-dominated  cases  (Rogers  et  al. 
20G3). 

3)  We  make  comparisons  at  only  one  location  (al  the 
location  of  buoy  45007  in  Fig.  1), 

4)  For  mode  I -data  comparisons,  we  use  a  location  near 
the  center  of  the  lake.  The  depth  is  165  m,  which  is 
relatively  deep  water  for  the  typical  wave  frequen¬ 
cies  in  the  lake.  Thus,  the  impact  of  finite  depth 
physics  is  limited. 


x  (km) 


FlO,  i.  Lake  Michigan,  with  depth  contours  (m)  and  NDBC 
instrument  locations  shown. 

4)  Challenge:  Describing  frequency  variation 

The  primary  challenge  with  quantitative  directional 
validation  of  a  long  time  series  is  that  a  different  set  of 
low-order  moments  exists  for  every  frequency  band. 
That  is  one  dimension.  Combine  that  with  the  time 
dimension,  and  the  validation  quickly  becomes  unman¬ 
ageable.  One  can  make  a  qualitative  comparison  by 
plotting  these  moments  as  a  function  of  time  and  fre¬ 
quency,  but  our  objective  is  to  make  quantitative  com¬ 
parisons.  Thus,  it  is  necessary  to  perform  some  kind  of 
integration  in  frequency  space.  Yet  we  cannot  throw 
out  the  frequency- wise  variation  of  these  moments  al¬ 
together,  since  (as  was  mentioned  in  section  1)  one 
objective  of  this  study  is  to  determine  whether  an  op¬ 
erationally  used  wave  model  adequately  reproduces  the 
directional  spreading  as  a  function  of  frequency  relative 
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to  the  spectral  peak.  Thus,  there  are  two  competing 
mot i valors:  I )  the  desire  lo  make  the  problem  more 
manageable  via  frequency- wise  integration  of  direc¬ 
tional  metrics  and  2}  the  desire  to  describe  the  fre¬ 
quency-wise  variation  in  directional  spreading. 

Our  approach  is  a  compromise  between  these  two 
motivators.  We  retain  frequency- wise  bins,  but  use 
fewer  bins  than  are  used  in  the  model  computational 
grid: 

1)  0. 5-0,8  times  the  relative  frequency  f/fp  ("low  fre¬ 
quencies"), 

2)  Q£f/fp-l2f/fft  {“frequencies  at  and  near  the  peak"). 

3 )  1 ,  2f/fp- 2 i }f/fp  ( ' *  f  req  ue  n  c  \  e  s  a  be  >  v e  t  h  e  pe  a  k  " ) .  a  n  d 
4}  2.()f/fp-3A)f/fp  (“highest  frequencies"), 

5)  CHALLENGE:  DEFINING  Jill  PEAK  FREQUENCY 

To  quantify  the  variation  of  the  directional  spreading 
as  a  function  of  relative  frequency,  ii  is  obviously  nec¬ 
essary  to  define  the  peak  frequency.  Though  this  may 
sound  simple,  it  is  subject  to  problems,  since  even  in  a 
region  like  Lake  Michigan,  with  its  typically  simple  sea 
stales,  peak  frequency  can  be  a  rather  unstable  quan¬ 
tity,  with  significant  model-data  mismatch  being  not 
uncommon.  Obviously,  it  is  very  problematic  to  com¬ 
pare  model  directional  spreading  as  a  function  of  mod¬ 
eled  relative  frequency  to  observed  directional  spreading 
as  a  function  of  observed  relative  frequency  in  cases 
where  the  modeled  and  observed  peak  frequencies  are 
very  dissimilar.  Model  predictions  of  mean  period  tend 
to  be  more  reliable  and  much  more  stable.  To  address 
this,  we  use  a  “synthetic  peak  period."  which  is  a  simple 
function  of  the  mean  period.  The  relation  is  determined 
using  a  simple  linear  regression  of  the  two  metrics  for 
the  time  period  of  the  hindcasl  described  in  section  5 
below.  The  mean  period  is  calculated  over  the  fre¬ 
quency  range  of  0.07-0.4  Hz.  For  the  modeled  values, 
the  result  of  the  regression  is 

Tp  =  l.2I42rmc.m  -  0,7126. 

For  the  buoy,  the  regression  is 

l),  =  1.23257^  -  0.7051. 

In  subsequent  discussions,  Tft  and  fp  refer  to  this  syn¬ 
thetic  peak  period  except  in  one  case  where  it  is  explic¬ 
itly  staled  that  the  “true"  peak  period  is  presented. 

6}  Challenge:  Avoiding  parametric  models 
AND  DAT  A  -  AD  A  FT  I VF.  METHODS 

It  was  mentioned  in  section  I  that  one  objective  is  to 
avoid  using  a  parametric  model  {e.g..  the  cos2li  model) 
or  a  data-adaptive  method  (e,g,,  MLE  and  MEM)  lo 


infer  directional  characteristics  from  the  buoy  data.  The 
solution  is  simply  to  use  only  variables  directly  ex¬ 
tracted  from  what  the  buoy  measures:  we  transform  I  he 
model  to  yield  quantities  analogous  to  what  the  buoy 
measures.  This  approach  has  been  taken  by  others:  for 
example,  Ardhuin  et  al.  (20113).  The  specific  calcula¬ 
tions  are  described  in  section  2h. 

h.  Definition  of  directional  metrics 

There  exists  a  separate  directional  distribution  tune- 
lion  for  each  frequency  component  that  can  he  decom¬ 
posed  into  a  Fourier  series: 

D\  f  tft  -  -j  -  +  ^  |«llcos(//<#l  4-  htf  sinudb |  j . 

(I) 

where 

a„(f)  ~  J  IM  f  (J)  cos/fO  do  and  (21 

hj  f)  -  J  D{  f  0)  sin#i0  tIU. 

J  n 

The  first  four  Fourier  coefficients  (u,.  M  can  be 

inferred  from  the  signals  measured  by  a  heave  -pitch 
roll  directional  buoy.  This  permits  only  an  approxima¬ 
tion  from  the  truncated  Fourier  series  ( Longue  i- 
Higgins  et  al,  1963;  Kuik  et  al.  1988),  which  is 

/}*(/,  0)  =  -|-+  2  [«!,(/) .cosfwfl)  +  hjf)  sim/Mb]  | 

(3) 

Unfortunately,  Eq,  (3)  has  limited  utility  for  describing 
D{  f  u),  since  it  is  only  accurate  if  the  unmeasured, 
higher-order  Fourier  components  are  very  small.  One 
possible  manifestation  of  this  inaccuracy  is  negative  val¬ 
ues  of  />*(£  f/).  Parametric  models  (such  as  (be  cm' 
form)  have  been  developed  to  yield  more  natural  (and 
thus  presumably  more  accurate)  representations  of 
D(f  0)  given  the  measured  low-order  moments,  but 
these  models  give  details  of  />(  j\  M)  that  are  not  actually 
determinable  from  buoy  motion.  Further,  al  least  one 
commonly  used  data-adaptive  method  -the  maximum 
likelihood  estimator— produces  a  D(  j\  0)  that  is  incon¬ 
sistent  with  the  original  cross-spectral  matrix  elements 
(Oltman-Shay  and  Guza  1984).  Kuik  el  al.  f  1988)  sug¬ 
gest  “model  free"  expressions  for  mean  wave  direction 
0(,  and  directional  width  <r„.  Kuik  el  ah  also  suggested 
two  higher-order  statistics  (skewness  and  k miosis)  that 
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we  do  not  use  herein.  All  four  statistics  are  expressible 
as  functions  of  the  four  Fourier  coefficients 
b\ (/)*  (h  (/)*  f>2(/)]-  Mean  wave  direction  is  given  as 
if)  =  arcum|/j,(/)/ff,(/)|.  Directional  width  is  quan¬ 
tified  as  the  '"circular  RMS  spreading*': 


(4) 


</)  =  V2[1 

where  »?,(/)  =  V«i(/)2  +  '?>, (_f)2. 

The  calculation  in  reverse  is  a,  =  m j  cos0(,  and  />,  — 
mx  sin% 

Real-time  and  historical  data  from  directional  Na¬ 
tional  Data  Buoy  Center  (NDBC)  buoys  include  esti¬ 
mates  of  the  low-order  moments  %  and  mx  (Steele  el  al. 
1985). 1  We  will  first  discuss  in  detail  the  calculation  of 
the  mean  direction  and  then  do  the  same  for  the  direc¬ 
tional  spreading. 

Calculation  of  Mean  wave  direction  and 

DIRECTIONAL  SPREADING 

In  the  literature,  the  mean  direction  is  the  most  com¬ 
monly  presented  directional  property  of  waves  (e.g.,  in 
maps  of  wave  heights  with  arrows  representing  mean 
direction).  Models  such  as  SWAN  (Booij  et  al.  1999} 
and  W  A  V  E  WATCH - 1 1 1  (WW3;  To!  man  1991,  2002) 
directly  calculate  actual  two-dimensional  spectra  £(/), 
ft)  and  out  put -averaged  %  for  frequencies  fx-f2  calcu¬ 
lated  as 


fy,  =  a  ret  an 


©■ 


(5) 


where 


B)  cq&B  df  dB/E  and 


,=  rr«/. 

J  n  J  r, 

f2* 

h  |  =  I  E(fMmBdfd(i/E 

Jo  Jr i 


(6> 


where  E  =  /,(’  E(  f.  6)dfdB.  (SWAN  and  WW3 

are  coded  to  output  fl,,  only  for/,  and  f2  equal  to  Hand 
respectively.) 

Now,  we  want  to  derive  %  based  on  u,.  hx  from  ihe 
buoy's  measurements  according  to  the  SWAN-WW3 
definition  in  (5),  We  use 


1  On  notation  used  elsewhere:  NDBC  uses  the  notation 
instead  of  “a/1  (used  by  Kuik  el  al.  and  herein)  and  “r/  instead 
of  "/«,*'  (used  by  Kuik  el  al.  and  herein).  Further,  ihe  NDBC 
definitions  of  ihe  Fourier  coefficients  (tf,.  />-»)  |as  used  in 

ihetr  literature  such  as  Steele  ct  al.  (1985)]  arc  dimensional, 
whereas  we  use  ihe  Kuik  convention  of  nundimensional  Fourier 
coefficients  (n,.  k?)*  The  notation  is  useful,  as  it  indi¬ 

cates  a  relation  to  (a,,  6,). 


where  /: 


a  i  =  f  a  i  (/)£(/)  df/E  and 
J  r, 

*,-f‘ 

J  f\ 

sT 


bdf)F{f)df!E, 
F(f)df : 


17) 


Note  that  if  wc  choose  fx  and  f2  as  values  close  t  of  say 
ft  =  0,9j fp  and  f2  —  \Afr  this  is  in  practice  very  similar 
to  the  "mean  wave  direction  corresponding  to  energy  of 
the  dominant  period'’  (MWD)  reported  by  NDBC  and 
to  the  "Dt"  reported  by  the  Coastal  Data  Information 
Program  (CD IP)  ("mean  direction  from  which  energy 
is  coming  at  the  peak  period").  The  use  of  a  broader 
band  of  frequencies  makes  the  metric  more  stable,  but 
increases  the  risk  that  two  distinct  wave  systems  could 
be  integrated  together. 

As  with  mean  direction,  the  directional  spreading 
trt,  —  o,i {/).  In  this  study,  we  use  a  weighted  mean  of  <r„ 
over  particular  frequency  ranges.  We  denote  i his  as  <th. 
To  be  consistent  with  our  calculation  of  mean  direction 
(5).  and  with  calculation  methods  of  SWAN  and  WW3, 
the  form  of  the  calculation  of  the  mean  directional 
spread  adopted  in  this  paper  is 


trn  =  ]2[  1  -Utr  +  b2)'2]  I1 
(WW3  uses/,  —  0  and  f2  =  *4. 


(8) 


With  the  model.  d}  and  /),  are  calculated  as  in  (6),  For 
buoy  measurements,  <7,  and  />,  are  calculated  as  in  (7), 

3,  An  idealized  case 

Rather  than  move  straight  to  the  hindeast  simulation, 
we  will  first  provide  an  idealized  application,  since  the 
idealized  application  is  used  as  a  point  of  discussion 
when  interpreting  the  hindeast  results. 

a.  Introduction  of  nonlinear  interaction 

cm npu tali* m  n methods 

One  limitation  of  the  dynamics  used  by  third- 
generation  (3G)  wave  models  is  the  highly  simplified 
D1A  (II asset mann  et  al.  1985)  used  to  compute  four- 
wave  nonlinear  interactions  in  both  models:  the  D1A 
uses  only  a  small  subset  of  the  possible  resonant  qua¬ 
druplets.  A  software  routine  based  on  the  Webb- 
Resio-Tracy  method  (WRT:  see  Resio  and  Perrie  1991 
and  references  therein)  has  been  implemented  in  the 
SWAN  model  by  G.  van  Vledder,  designated  as  "XnP 
in  the  user*s  manual.  In  contrast  to  the  DIA,  ill  is 
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Fig*  2,  Nondirecticmal  spectral  density  distributions  for  the  ide¬ 
alized  simulations  and  circular  rms  spreading  for  the  simulations 
with  S  =  The  vertical  lines  in  the  bottom  panel  indicate  the 
peak  frequencies  of  each  model. 


Flu.  3.  Directional  spectral  density  distributions  for  the  S  SnlA 
idealized  simulations:  (lop)  the  WRT  (Xnl)  and  (hoi tom)  the 
D1A  results. 


method  is  essentially  exact,  hut  relatively  time  consum¬ 
ing. 

h.  Simulation  descriptions 

An  example  application  of  this  WRT  subroutine  is 
shown  in  Fig.  2.  The  computation  is  with  a  “point 
model/'  implying  either  zero  propagation  or  infinite 
fetch:  (dCMtNldx)  +  (dC^N/fiy)  =  l).  First,  a  “spin up” 
simulation  was  run  using  all  three  deepwater  source 
terms  |D1A  for  nonlinear  interactions.  Tolman  and 
Chalikov  (1996)  for  wind  input  and  dissipation],  a  con¬ 
stant  wind  speed  of  Uni  —  18  ms  and  a  duration  of 
1  day*  The  resulting  spectrum2  was  used  to  initialize  two 
simulations  that  are  identical  except  that  one  uses  WRT 
and  the  other  PI  A.  These  two  simulations,  also  of 
1-day  duration,  include  only  nonlinear  interactions,  to 
lend  insight  regarding  the  effect  of  nonlinear  interac¬ 
tions  on  swell  as  it  leaves  its  source.  To  summarize,  the 
following  spectra  are  presented  here: 

1)  the  initial  conditions  for  the  other  two  simulations, 
which  is  the  final  condition  of  the  spin  up  simulation, 
and  which  includes  all  three  deepwater  source-sink 
terms.  S  =  Sin  +  5*  +  S„l4; 

2)  the  final  condition  of  a  simulation,  which  includes 
only  four-wave  interactions,  S  —  Sf,]4,  calculated  us¬ 
ing  the  WRT  routine:  and 


The  shghl  difference  between  the  two  initial  conditions  is  due 
to  the  difference  in  frequency  resolution  for  the  two  simulations: 
with  the  D1A  method,  a  logarithmic  distribution  with  f,  1.1/)  , 
is  preferred  and  with  the  WRT  method,  a  higher  resolution  is 
preferred:  we  use  f]  =  L 07/  ,  (4S  frequency  bins  from  (MW  18  to 
1,0  Hz). 


3)  the  final  condition  of  a  simulation,  which  includes 
only  four-wave  interaction,  S  =  5nW,  calculated  us¬ 
ing  the  D1A  routine. 

Note  that  since  this  model  does  not  include  propaga¬ 
tion,  dispersion  of  the  swell  is  not  represented.  The 
effects  of  dispersion  could  be  significant  within  I  day, 
depending  on  the  size  of  the  storm:  this  dispersion 
would  be  expected  to  reduce  wave  steepness  and  there¬ 
fore  nonlinear  interactions*  The  differences  seen  here 
between  the  PI  A  and  WRT  models  are  qualitatively 
consistent  with  compulations  of  the  nonlinear  source 
term  by  Hasselmann  el  al  ( 1985;  e.g*.  see  their  Fig.  7), 

c.  D  isci tssi on  c if  rest  i  its 

Figure  2  shows  the  nondireetional  spectral  density  of 
the  three  spectra  (top  panel)  and  the  directional 
spreading  of  the  second  and  third  spectra  (bottom 
panel).  Since  only  three  spectra  arc  being  presented 
with  no  time  dimension,  it  is  not  necessary  to  integrate 
in  frequency  space,  and  the  actual  variation  with  fre¬ 
quency  at  the  model  resolution  is  shown*  Skewness  and 
kurtosis  for  the  second  and  third  spectra  were  also  com¬ 
pared,  but  the  comparisons  were  not  noteworthy  and 
are  not  presented  here.  The  two-dimensional  spectral 
density  distributions  for  the  simulations  with  S-  -  Snu 
are  shown  in  Fig.  3,  In  this  figure,  both  spectra  have 
been  normalized  by  1.19  nr  Hz  1  \  which  is  the 

maximum  of  the  third  spectrum.  Thus,  the  contours  are 
labeled  relative  to  the  peak  of  the  larger  spectrum.  The 
fo  1 1  o w i  ng  o  bse  r v at  i  ons  ea  n  be  ma  d  c . 

*  Though  it  is  not  directly  related  to  the  subject  matter 
of  this  study,  the  effect  of  the  inaccuracy  of  the  PI  A 
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on  frequency  downshifting  is  seen  clearly  in  the  spec¬ 
tral  density  plot. 

*  The  inaccuracy  of  the  D1A  is  leading  to  slightly  too 
narrow  a  directional  distribution  at  the  low  frequen¬ 
cies,  slightly  too  broad  a  distribution  near  the  peak, 
and  clearly  too  broad  spectra  above  the  peak,  most 
noticeable  beyond  0.1  Hz. 

*  The  directional  spectrum  plot  (Fig.  3)  gives  an  imme¬ 
diate  visual  impression  that  directional  spreading  is 
much  greater  with  the  DIA  model;  this  is  reflected 
quantitatively  in  Fig,  2.  However,  the  higher  direc¬ 
tional  spreading  is  really  apparent  only  in  the  lowest 
energy  contour  (2%  of  peak),  at  the  higher  frequen¬ 
cies, 

4,  A  hindcast  validation 

a.  Simulation  description 

The  grid  domain  is  shown  in  Fig.  I.  The  following 

settings/features  were  identical  to  those  of  Rogers  el  al, 

(2(K)3). 

*  Cartesian  coordinates  were  used,  with  grid  spacing  of 
2  km, 

*  The  lake  bathymetry  is  provided  by  the  National 
Oceanic  and  Atmospheric  Administration  (NOAA)/ 
Great  Lakes  Environmental  Research  Laboratory, 

*  The  directional  resolution  is  10°. 

*  The  wind  field  is  created  using  wind  observations 
from  the  two  open-water  buoys  in  Lake  Michigan 
{45002  and  45007),  adjusted  to  LO-m  elevation,  with 
linear  interpolation  in  the  latitudes,  and  no  variation 
in  longitude, 

*  Default  param eterizations  for  Sm,  Sdsi  and  SnU  are 
used,  except  that  the  power  on  the  relative  wavenum¬ 
ber  [denoted  //  in  Rogers  ct  ah  (2003)]  is  set  to  2.0. 
( The  default  para  me  terizat  ions  in  SWAN  are  that  of 
W  AM,  cycle  3,  sometimes  referred  to  in  the  literature 
as  "W  A  M3  physics,”) 

The  following  settings/features  are  different  from 

those  of  Rogers  et  al,  (2003). 

*  Season  hindcast  covers  tXMJO  UTC  1  September-0500 
UTC  14  November  2002, 

*  The  frequency  grid  is  logarithmic,  with  29  frequencies 
from  0.07  to  L0  Hz. 

*  Since  l  September  2002  was  relatively  calm,  only  a 
very  short  ''ramp”  time  was  needed  (6  h),  so  the  com¬ 
parisons  to  the  data  start  at  0600  UTC  1  September, 

*  A  time  step  of  6  min  is  used. 

*  The  version  of  SWAN  used  is  40,41  A  (Booij  et  al. 
2005). 
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The  physical  pa  ram  eterizations  used  are  not  tuned 
for  this  simulation  or  for  this  area:  rather,  they  are  the 
same  as  those  used  in  the  SWAN  forecasting  systems 
run  at  NRL  for  other  areas. 

h.  Simulation  results 

The  primary  focus  of  the  study  is  the  accuracy  of  the 
predictions  of  directional  spreading  in  the  hindcast. 
However,  results  of  any  validation  of  directional 
spreading  will  be  much  more  meaningful  if  it  is  first 
shown  that  the  nondirectionai  spectra  and  the  mean 
direction  are  well  predicted.  Thus,  we  present  results 
other  than  the  directional  spreading  before  making  the 
comparisons  of  directional  spreading, 

1)  Results:  Nondirectionai.  spectra  and 

MEAN  WAVE  DIRECTIONS 

To  provide  a  sense  of  the  length  of  the  simulation 
and  how  many  events  are  being  verified,  a  time  series  of 
zero-moment  wave  height  //,„<„  at  buoy  45007,  is  shown 
in  Fig,  4,  These  wave  heights  are  also  compared  to  data 
in  scat  ter  plot  form,  along  with  the  mean  period,  the 
mean-mean  wave  direction,  and  the  true  peak  period, 
in  Fig.  5?  Statistics  associated  with  the  comparison  are 
indicated  in  the  plots.  The  wave  height  and  mean  pe¬ 
riod  arc  for  the  frequency  range  0,07-0,4  Hz  (essen¬ 
tially  the  entire  spectrum).  The  mean-mean  wave  di¬ 
rection  is  the  mean  wave  direction  integrated  over  0.8 
fp-l2fp  using  (5),  so  it  is  a  stable  metric  of  the  mean 
direction  near  the  peak  frequency.  By  the  standards  of 
a  wave  model  that  uses  only  wind  forcing,  the  agree¬ 
ment  is  very  good  for  all  four  metrics.  The  good  pre¬ 
diction  of  the  wave  height  and  mean  period  suggests 
that  the  nondirectionai  wave  spectra  F{f)  are  fairly 
well  predicted.  This  provides  confidence  that  the  hind¬ 
cast  is  suitable  for  detailed  study  of  the  accuracy  of  the 
prediction  of  the  directional  spreading.  Some  bias  is 
evident  in  the  prediction  of  the  mean  period,  indicating 
a  problem  with  overestimation  of  energy  below  the 
peak  or  underestimation  of  energy  above  the  peak,  or 
both. 

Even  with  excellent  agreement  in  the  three  nondi- 
reetional  parameters,  there  can  still  be  problems  with 
the  frequency  width  that  are  not  revealed.  Thus,  we 
present  in  Fig.  6  time -col  located  scat  ter  plot  compari¬ 
sons  of  the  energy  level  in  the  four  frequency  bands 
described  in  section  2,  quantifying  the  bias  and  random 


■  In  the  ease  of  the  peak  period,  the  density  of  occurrence  is 
plotted  rather  than  individual  points,  since  discrete  peak  period 
values  tend  to  overlay  each  other. 
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error  in  each  hand*  The  "partial  wave  height'*  pre¬ 
sented  in  Fig.  6  is  calculated  from  ilie  variance  (i.e.* 
energy)  of  the  wave  spectrum  over  a  frequency  range 
defined  by  lower  and  upper  hounds  /,  and  /*: 

--  4\/Vrciui  and  Ppiiriini  =  F(f)4f i  the  "partial  vari- 
a  nee,"  (The  fictitious  quantity  is  used  rather 

than  variance,  since  wave  height  is  more  intuitively  un¬ 
derstood.)  A  time-averaged  non  direct  ion  a  I  spectrum 
F(  f/fjJ  is  shown  in  Fig,  7*  This  is  created  hy  using  24/// 
bins  instead  of  4.  Since  it  is  time  averaged*  it  quantifies 
bias  only*  The  robust  feature  in  Figs,  6  and  7  is  an 
overestimution  of  energy  below  ihe  spectral  peak,  thus 
explaining  the  modest  positive  bias  in  the  mean  period* 
Interestingly*  the  model-data  comparison  here  is  quali¬ 
tatively  very  similar  to  the  F(f)  comparison  for  ihe 
idealized  case  (Fig*  2*  top  panel)*  This  similarity  sug¬ 
gests  that  the  overprediction  of  low-frequency  energy 
in  Fig*  7  is  at  least  partially  attributable  to  inaccuracy 
associated  with  the  DIA,  The  problem  can  be  compen¬ 
sated  for  by  reducing  the  weighting  on  the  relative 
wavenumber  in  the  Komen  et  al.  ( I9S4)  ,Vds  formulation 
used  by  SWAN,  but  this  would  only  shill  the  positive 
bias  to  the  frequencies  above  the  peak:  this  has  been 
verified  by  repealing  the  hindeast  with  a  weighting  of 
n  =  1*5  instead  of  n  =  2*0*  jin  SWAN,  the  default 
setting  is  n  ~  1.0,  but  this  setting  consistently  leads  to 
underprediction  of  the  mean  period  in  cases  of  wind 
speeds  up  to  21  m  s  J;  for  more  detailed  discussion  of 
this  tuning  parameter  in  SWAN,  see  Rogers  et  al* 
(2003)*] 


In  summary,  low-frequency  energy  is  over  predicted 
by  the  model  in  the  hindeast  (bias  -  9  cm.  r  0*04), 
and  ibis  should  be  considered  when  evaluating  direc¬ 
tional  spreading  in  this  frequency  range* 

2)  Results:  Directional  spreading 

{ i )  Sc  ratterp  h  >  t  cot  t  ip  a  rixot  tx 

The  scatlerplot  comparisons  of  the  mean  directional 
spreading  <rtt  arc  made  in  Figs,  ft  a  and  8b.  Figure  ft  a  is  a 
simple  scatlerplot  comparison  of  <rfll:  statistics  associ¬ 
ated  with  the  comparison  are  shown  in  each  plot.  In  Fig* 
8b*  the  horizontal  axis  is  the  buoy  partial  wave  height 
for  the  indicated  frequency  range*  and  the  vertical  axis 
is  the  misfit  in  the  mean  directional  spreading,  <r„,K.  - 
trHuhs.  There  are  fewer  points  in  the  highest -frequency 
comparisons  (2/  to  3/)  because  the  highest  frequency 
in  the  directional  buoy  data  is  0.35  Hz:  thus,  there  are 
often  no  data  available  in  this  frequency  range,  depend¬ 
ing  on  the  value  of  /,  In  alt  plots  of  the  directional 
spreading,  "weak  signal"  data  points  are  not  included, 
being  defined  as  collocated  values  for  which  either  ihe 
buoy  or  modeled  total  wave  heights  (Fig*  5)  are  less 
than  0.5  m. 

We  make  the  following  observations* 

•  Low  free p  t  end  ex  ( ( 1,  6  /-( I .  ft/ ) :  S  W  A  N  u  n  de  r  p  re  diets 
spreading  (bias  of  -24  )  and  ihe  re  is  much  scalier 
\erms  ~  27  and  the  scatter  index  (SI)  -  (1.52 1,  Note 
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0  of  buoy  (  °) 


mean  period  of  buoy  (s) 


Flo.  5.  Scaitcrplot  comparison  of  wave  height,  mean  period,  true  peak  period,  and  mean-mean  wave  direction 
for  the  hindcasi.  The  statistics  listed  are  correlation  coefficient  (r).  scatter  index  (SI),  sid  dev  of  error  (ksD).  rms 
error  and  mean  error  (bias).  In  the  plot  of  true  peak  period,  points  are  not  plotted'  instead,  the  density  of 

points  is  indicated  by  the  si^e  of  the  circles. 


that  in  this  case,  even  the  mean  value  tor  the  "ground 
truth"  is  not  reliable:  this  is  discussed  in  detail  in 
section  6. 

*  Freqt ter. i ctes  # 1 car  p eak  ( 0 . 8 f  - \2fr):  Ra ndo m  e rro r  i s 
smaller  (eRMS  -  6.5  ),  but  still  not  as  good  as  it  is  for 
the  other  metrics  (wave  height,  etc.).  There  is  not  a 
discernable  systematic  error  (bias  -  1.2"),  The  agree¬ 
ment  is  especially  good  for  moderate  and  large  wave 
heights  (Fig.  hh).  Note  that  the  buoy  data  are  more 
reliable  for  these  moderate  and  large  wave  heights 
(Anctil  et  al.  L993). 

*  Frequencies  above  the  peak  ( 1 2fp-2fp):  SWAN  does 
not  do  a  very  good  job  of  following  the  observations 
(predicted  spreading  varies  much  less  than  the  ob¬ 


served  spreading,  r  —  0,44),  but  error  tends  to  he  low4 
(fHMS  -  5.6°),  and  there  is  no  significant  systematic 
error  (bias  =  2°). 

*  Highest  frequencies  (2fp-3f„):  Like  the  prior  fre¬ 
quency  range.  SWAN  does  not  do  a  very  good  job  of 
following  the  observations  (r  =  0,44):  predicted 
spreading  is  consistently  close  to  40".  However,  again 
the  error  lends  [a  be  low  (eRMS  =  4.5"),  since  the 


4  Judgment  of  what  constitutes  “low"  rtm  error  in  directional 
spreading  is  necessarily  subjective:  it  can  be  compared  with  ex¬ 
pected  measurement  uncertainty,  for  substantiation — why  this 
level  of  random  error  might  be  considered  “low” — sec  section  6, 
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Fir.,  6,  Comparison  of  “partial  wave  height”  {m}.  describing  the  energy  level  in  the  four  frequency  bands 

described  tn  section  2. 


observations,  though  they  show  more  variation,  are 
also  clustered  near  40°, 

(ii)  Time-averaged  comparisons 

To  perform  time  averaging,  the  hindeast  and  ob¬ 
served  directional  spreading  are  calculated  over  smaller 

frequency  bins  of  0Afp  (so  the  bins  are  0.5/,,,  0.6/ . , 

2.7  /,,  2.8/,).  To  enhance  stability,  the  integration  to 
calculate  tr„  is  performed  over  a  ±tKl  fp  range,  so  points 
are  used  more  than  once,  similar  to  a  moving  average 
comparison.  A  simple  time  averaging  is  used  (Le„  the 
values  are  not  weighted).  The  resulting  distributions 
are  shown  in  Fig.  9,  along  with  the  empirical  parametric 
model  of  Donelan  et  at,  ( 1985),  which  was  extended  by 
Banner  (1990)  [see  also  Young  (1999),  Eq.  (5.66) |. 

At  the  lower  frequencies,  directional  spreading  of  the 


buoy  is  approximately  60%  higher  than  that  of  either 
the  parametric  model  or  the  numerical  model.  At  the 
highest  frequencies,  the  directional  spreading  of  the 
parametric  model  is  approximately  25%  higher  than 
that  of  the  buoy  and  the  numerical  model. 

5.  Summary  uf  results  regarding  bias  in 
directional  spreading 

Our  reading  of  the  literature — specifically,  Khande- 
kar  el  al.  (1994),  Forristall  and  Ewans  (1998),  Furr ist all 
and  Greenwood  (1998).  Cardone  and  Resit)  (1998), 
Wyatt  et  al.  (2003)— gave  us  the  impression  that  third- 
generation  wave  models  such  as  WAM  have  a  fairly 
consistent  tendency  to  over  predict  directional  spread¬ 
ing.  In  Forristall  and  Greenwood  (1998)  (and  see  also 
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March  2<K>7  - 


Fie.  7.  ( Lop)  Time-averaged  non  directional  spectrum  F(f/fp). 
(bottom)  Number  of  time  records  used  to  create  a  time  average  at 
each  of  the  24  f/f{t  bins. 

Cardone  and  Resio  1998),  the  problem  is  quite  reason¬ 
ably  attributed  to  inaccuracy  associated  with  the  DIA 
approximation  of  nonlinear  interactions  used  in  third- 
generation  wave  models* 

In  this  section  we  consider  the  results  of  the  hindcast 
together  with  those  of  the  idealized  simulations  and 
contrast  both  to  our  prior  expectations*  For  the  ideal¬ 
ized  simulations,  a  model  with  exact  calculations  of 
nonlinear  interactions  is  taken  as  "ground  truth"  while 
for  the  hindcast,  buoy  data  are  taken  as  ground  truth. 

The  comparison 

Low  frequencies:  Our  prior  expectation  was  that  the 
model  directional  distribution  would  be  too  broad* 
but  in  the  long  hindcast,  the  model  directional 
spreading  ts  narrow  relative  to  the  ground  truth.  In 
the  idealized  case  the  "model"  is  close  to  the 
ground  truth  (only  slightly  too  narrow). 

Near  the  peak:  Our  prior  expectation  was  that  the 
model  directional  distribution  would  be  too  broad, 
but  in  both  the  idealized  case  and  the  long  hind- 
cast*  the  average  model  directional  spreading  is 
quite  close  to  that  of  the  ground  truth* 

High  frequencies:  Our  prior  expectation  was  that  the 
model  directional  distribution  would  be  too  broad. 
The  idealized  simulation  supports  this,  but  in  the 
long  hindcast,  the  average  model  directional 
spreading  is  quite  close  to  that  of  the  ground  truth. 

As  a  sort  of  a  disclaimer,  we  refer  above  to  other 
third-generation  wave  models  used  in  prior  studies*  We 
do  not  imply  that,  had  we  used  another  3G  wave  model 
tn  our  hindcast,  our  results  would  be  the  same.  We  are 
contrasting  our  results  with  our  prior  expectations*  Just 


as  biases  in  predictions  of  nondirectional  moments 
(e.g.*  wave  height*  mean  period)  are  obviously  sensitive 
to  a  model's  source-sink  formulations  (wind  input,  dis¬ 
sipation*  etc.)*  biases  in  predictions  of  directional  width 
rrH(  f )  should  be  expected  to  be  sensitive  to  these  for¬ 
mulations. 

6*  Discussion 

a .  Accuracy  of  mean  direction  in  turning  winds 

Comparisons  of  mean  wave  directions  for  this  hind- 
cast  simulation  show  rather  good  accuracy  overall.  The 
response  of  a  3G  wave  model  to  rapidly  turning  winds 
is  a  concern.  We  do  not  specifically  address  this  prob¬ 
lem  here*  but  we  do  not  mean  to  imply  that  it  is  not  an 
area  in  w'hich  the  models  may  bear  significant  improve¬ 
ment. 

b.  Fetch  geometry 

The  influence  of  the  fetch  geometry  is  represented 
within  the  formulation  of  third-generation  wave  models 
such  as  SWAN,  Thus*  it  is  presumed  that  the  observed 
and  modeled  spectra  are  both  influenced  by  the  geom¬ 
etry  of  Lake  Michigan*  Ataktiirk  and  K  a  Isa  ros  (1999.  p. 
643),  cite  a  large  reduction  in  wave  energy — modeled 
and  observed — associated  with  the  narrow  width  of 
Lake  Washington*  Since  it  is  a  relatively  large  lake,  this 
effect  is  expected  to  be  much  less  pronounced  in  Lake 
Michigan,  perhaps  more  comparable  to  the  Lake  On¬ 
tario  observations  cited  in  Ataktiirk  and  Kalsaros 
(1999)*  which  are  from  Donelan  el  al.  (1985)* 

c*  The  challenge  of  mixed  seas  and  swells 

In  the  case  of  mixed  seas  and  swells*  the  challenge  of 
directional  validation  is  much  greater*  The  location  of 
our  hindcast  was  deliberately  chosen  to  avoid  this  ad¬ 
ditional  complexity.  Frequency -wise  integration  intro¬ 
duces  the  risk  of  mixing  multiple  components  (e.g.,  seas 
and  swells).  A  frequency- integrated  metric  (e.g.*  mean 
direction  or  directional  spreading)  that  includes  mul¬ 
tiple  wave  systems  is  of  dubious  value*  Under  such  cir¬ 
cumstances,  in  order  to  present  the  type  of  comparison 
made  here  (e.g.,  in  Figs.  6  and  7),  the  windsea  compo¬ 
nent  must  be  identified  and  separated  from  the  swell 
components* 

Methods  for  separating  individual  sea  and  swreil  com¬ 
ponents  in  a  wave  spectrum  exist  in  the  I  i  to  rain  re  (e.g.* 
Beal  1991;  Ceding  1992;  Komen  el  al.  1994;  Hanson 
and  Phillips  2001)*  Thus*  it  is  theoretically  possible  to 
compare  measured  and  observed  wave  spectra  in  a 
component- wise  fashion  ( Hanson  and  Jensen  2004). 
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Flu,  8,  (a)  ScaUurpksi  comparisons  {hiudcast  vs  observation)  of  mean  directional  spreading  (  )  over  four  fre¬ 
quency  ranges,  (h)  Scaltcrploi  comparisons  (iiindcast  vs  observation)  for  lorn  frequency  ranges.  The  horizontal  axis 
is  the  buoy  partial  wave  height  for  the  indicated  frequency  range,  and  the  vertical  axis  is  the  misfit  in  mean 
directional  spreading  ( hindcast  -  observed). 


Unfortunately,  due  to  model  limitations,  it  is  not  un¬ 
common  to  have  a  swell  system  ilia!  exists  in  observa¬ 
tions  hut  not  in  the  model  spectra,  or  vice  versa.  In  this 
ease,  validation  of  the  directional  spreading  is  obviously 
not  possible. 

Based  on  our  experiences,  we  expect  that  a  valida¬ 
tion  such  as  was  performed  here  would  be  difficult  for 
an  exposed  coastline,  with  frequent  mixed  sea-swell 
conditions.  In  such  a  case,  some  compromise  is  prob¬ 
ably  necessary.  By  way  of  summary,  two  possible  com¬ 
promises  are  to  either 

I }  consider  a  shorter  time  period,  so  that  qualitative 
comparisons  can  be  made,  for  example  by  graphing 
trtl  —  frM( /,  0  and  E  —  £(/,  r),  or 


2)  separate  the  windsea  from  swell  components,  and 
validate  the  directional  spreading  of  the  windsea 
component. 

Of  course,  it  is  possible  to  evaluate  the  directional 
spreading  of  swell  components  also.  This  was  done  tor 
observational  data  by  Ewans  (2(M)I).  However,  if  the 
objective  is  to  evaluate  generation-stage  source  sink 
terms,  study  of  the  directional  spreading  ol  the  v\  indsca 
seems  to  be  the  most  direct  approach, 

(I.  Sensitivity  of  results  to  method  of  calculation  of 
peak  frequency 

It  is  well  accepted  (e,g..  Donelan  et  a l,  19S5:  Banner 
1990)  that  a  typical  windsea  directional  distribution  will 
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tend  to  he  most  narrow  somewhere  near  the  peak  fre¬ 
quency*  The  most  natural  way  to  analyze  this  behavior 
(or  lack,  thereof)  in  a  model  and  collocated  observa¬ 
tions  is  via  normalization  of  the  results  by  peak  fre- 


FlO.  V,  Comparison  or  time-averaged  results  (model  and  obser¬ 
vation)  with  the  parametric  model  of  Don  elan  et  a).  (1985)  and 
Banner  (1990).  The  hot  tom  panel  shows  the  number  of  time 
records  used  to  create  a  lime  average  for  the  model  and  buoy  data 
at  each  of  the  17  f/ff,  bins. 


quency.  as  is  done  here.  However,  this  does  introduce 
some  ambiguity-  There  are  multiple  options  regarding 
which  peak  frequency  to  use  (model  peak,  buoy  peak, 
or  some  combination),  and  one  can  expect  that  results 
will  demonstrate  some  sensitivity  to  this  choice  unless 
fairly  broad  frequency  bands  (e.g,.  the  lour  bands  used 
in  Figs.  6  and  X)  are  used.  With  small  frequency  bands 
very  near  the  spectral  peak  (sav.  0*95-1.05 
one  should  not  be  surprised  if  even  the  sign  of  the  bias 
is  not  robust.  Thus,  it  is  critical  to  use  consistent  meth¬ 
ods  when  intercomparing  statistics  from  different  hind- 
casts. 

e.  Consideration  of  uncertainty  in  observations 

It  is  useful  to  put  the  model-data  misfit  in  the  context 
of  measurement  error.  Quoting  Voorips  et  aL  (1997). 
"Observation  errors  consist  of  instrumental  errors,  er¬ 
rors  due  to  the  random  variability  of  the  spectrum  and 
limited  sample  time,  and  representation  errors.  Repre¬ 
sentation  error  is  the  difference  between  what  the  buoy 
actually  measures  . . .  .  and  its  model  equivalent/'  Sta¬ 
tistical  error — determined  by  degrees  of  freedom — has 
been  successfully  incorporated  into  a  validation  of  non- 
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directional  spectra  by  Alves  el  al.  (2002),  Unfortu¬ 
nately,  no  such  method  exists  that  accounts  for  the  tola! 
error,  since  this  is  not  defined,  at  least  not  in  the  context 
of  spectral  density  from  3-m  discus  buoys, 

Kuik  el  al.  (1988)  estimate  the  confidence  limits  on 
directional  moments  based  on  statistical  uncertainty:  in 
terms  of  rms  error,  they  are  5  I  O'  for  mean  direction, 
]0%-J5%  for  directional  width,  30%-S0%  for  skew¬ 
ness,  and  25%-l(X)%  for  kurlosis  (see  also  Anciil  el  al. 
1993),  So,  0.10-0.15  can  be  compared  with  the  scatter 
index  values  given  in  Fig.  8a, 

Measurement  of  long,  low-amplitude  waves  by  buoys 
is  problematic  due  to  the  weak  acceleration  or  slope 
signal  against  the  background  noise.  However,  we  do 
not  know  of  any  arguments  that  this  may  manifest  itself 
as  bias  in  the  spectral  density.  In  contrast,  it  has  been 
shown  by  Kuik  el  al.  ( 1988)  that  low  levels  of  noise  in 
the  surface  elevation  or  slope  will  cause  a  positive  bias 
in  the  directional  spread  (see  also  discussions  in 
O'Reilly  et  al.  1996),  In  the  case  of  the  low  frequencies, 
the  directional  spreading  of  the  buoy  used  as  ground 
truth  in  the  hindcast  herein  is  almost  certainly  loo  high. 
This  behavior  is  consistent  with  the  overestimating  of 
the  directional  spreading  of  swells  by  NDBC  3-m  dis¬ 
cuss  buoys,  as  reported  by  O'Reilly  et  al.  (19%)  (a  b 
bias,  with  the  metric  integrated  over  0.06-0.14  Hz). 

To  summarize,  in  the  model-buoy  comparisons 
herein,  bias  below  the  spectral  peak  is  evident  both  in 
the  spectral  density  F(f)  and  the  directional  spread 
</„(/).  Hias  in  <rfj(/)  below  the  peak  cannot  he  defini¬ 
tively  attributed  to  the  model,  whereas  bias  in  F(f) 
below  the  peak  is  credibly  attributed  to  the  model. 

f  The  high-frequency  cutoff  in  the  idealized 

simulations 

The  WAVEWATCH-m  and  WAM  models  employ 
a  diagnostic  tail  above  a  prescribed  frequency  relative 
to  the  spectral  peak  frequency.  Banner  and  Young 
(1994)  point  out  that  removal  of  this  tail  has  dramatic 
consequences  on  all  quantities  derived  from  the  wave 
spectrum.  SWAN,  however  does  not  employ  a  self- 
adjusting  high-frequency  cutoff,  and  in  the  idealized 
simulations  herein,  we  use  a  high  frequency  cutoff  fixed 
at  1.0  Hz,  which  is  in  fact  higher  than  the  fixed  cutoff 
frequency  used  by  Banner  and  Young  (1994).  Thus* 
there  is  less  concern  about  the  effect  of  the  parametric 
tail  on  the  results  presented.  However,  these  same 
simulations  were  performed  with  the  WAVEWA1  CH¬ 
IU  model  (not  shown  herein),  and  the  qualitative  im¬ 
pact  of  the  nonlinear  solver,  WRT  versus  DIA.  in  this 
simulation  is  similar  regardless  of  which  model  is  used 
as  the  platform. 


g.  The  impact  of  the  nonlinear  solver  it i  a  hind  cast 

Though  we  apply  the  WRT  nonlinear  solver  in  an 
idealized  scenario,  it  would  be  possible  to  apply  it  in  a 
shortened  version  of  our  Lake  Michigan  hindcast  to 
specifically  study  the  impact  of  the  inaccuracy  of  the 
DIA.  Presuming  that  the  WRT- based  model  would 
have  narrower  spreading  in  high  frequencies  (vs.  the 
DIA-hased  model),  we  can  reasonably  expect  that  the 
WRT-based  method  would  underpredict  the  directional 
spreading  in  the  high  frequencies  in  this  hindcast.  Ob¬ 
viously.  this  indicates  a  situation  in  which  a  model  (the 
DIA -based  model)  is  correct  for  the  wrong  reasons. 
Indeed,  if  this  pattern  is  systematic  (observed  in  other 
hindcast s).  it  would  justify  a  retiming  of  the  directional 
spreading  of  the  wind  input  term  in  SWAN  to  create 
broader  directional  spreading  in  the  high  frequencies 
with  the  WRT-based  model.  We  believe  that  in  any 
case,  a  move  to  more  accurate  calculations  of  nonlinear 
interactions  will  necessitate  retiming  of  the  other  two 
deepwater  source-sink  terms. 

In  fact,  a  hindcast  with  a  WRT-based  model  has  been 
performed  recently  (F.  Ardhuin  2(KJ5,  personal  commu¬ 
nication),  which  suggests  that  DI  A  does  lead  to  broader 
directional  spectra  in  the  higher  frequencies,  compared 
to  a  model  with  exact  nonlinear  computations.  With  the 
WAM3  physics  (also  employed  in  the  present  study), 
the  DIA-based  model  tends  to  be  too  broad  overall  and 
the  WRT-based  model  lends  to  be  too  narrow,  in  the 
higher  frequencies,  compared  to  observations.  Further, 
an  interesting  question  is  raised  by  the  authors  of  that 
study:  whether  modeled  directional  spreading  is  sensi¬ 
tive  to  directional  spreading  of  the  ,SMI  term,  or  con¬ 
trolled  solely  by  the  Snl4  term.  Further,  our  experience 
is  that  modeled  directional  spreading  is  sensitive  to  the 
.VtK  term:  we  have  verified  that  out  a  priori  choice  of  the 
weighting  on  the  relative  wavenumber  in  Sfc  (section 
4a)  does  affect  the  bias  statistics  presented  in  Fig.  8. 
Arguments  exist  in  the  literature  that  high-frequency 
directional  spreading  in  nature  is  controlled  solely  by 
the  Snl 4  term  (Young  and  Van  Vledder  1993:  Banner 
and  Young  1994:  Young  el  al,  1995). 

h.  Potential  subsequent  work 

A  subsequent,  more  ambitious  directional  validation 
exercise  might  he  of  longer  duration,  with  multiple  ob¬ 
servational  points,  and  might  also  consider  higher- 
order  moments:  skewness  and  kurlosis. 

7*  Conclusions 

In  an  enclosed  basin  such  as  Lake  Michigan,  it  is 
demonstrated  herein  that  il  is  possible  to  quantitatively 
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validate  directional  characteristics — mean  direction 
and  directional  spreading — of  a  relatively  long  wave 
model  hind  cast.  Further,  huoy  observations  can  be  used 
in  such  a  validation  without  applying  a  parametric 
model  (such  as  the  cos2s  model)  or  a  data-adaptive 
method  (such  as  the  MLE  or  MEM)  to  the  observa¬ 
tions.  Populations  of  model-observation  pairs  such  as 
(he  scatterplot  comparisons  herein  are  readily  con¬ 
densed  to  statistics  such  as  root-mean -square  error, 
bias,  and  standard  deviation  of  error,  so  it  is  feasible  to 
present  directional  validations  for  multiple  locations 
within  limited  space,  such  as  a  journal  article. 

Due  to  the  considerable  added  complexity  associated 
with  mixed  sea-swell  conditions,  it  is  not  as  straightfor¬ 
ward  to  perform  a  validation  in  this  manner  on  an  ex¬ 
posed  coastline.  Some  sea-swell  separation  algorithm 
would  need  to  be  applied. 

In  addition  to  the  validation  of  the  long  bindcast,  a 
pair  of  idealized  simulations  is  presented  herein  to  aid 
tn  the  interpretation  of  the  results.  The  two  idealized 
simulations  differ  in  their  methods  of  calculating  the 
nonlinear  interaction  [essentially  an  exact  method 
(WRT)  versus  an  operational  method  ( D I A)}_  Consid¬ 
ering  both  the  hindeast  validation  and  the  idealized 
simulations,  we  find  the  following  about  the  bias  char¬ 
acteristics. 

•  At  frequencies  below  the  spectral  peak,  in  both  the 
idealized  case  and  the  long  hindeast,  Lhe  model  di¬ 
rectional  spreading  is  narrow  relative  to  the  ground 
truth.  In  the  case  of  the  hindeast,  however,  the  signal - 
to- noise  ratio  of  the  ground  truth  may  be  poor  at 
these  low  frequencies.  In  the  example  of  the  idealized 
case,  the  bias  is  very  slight.  In  both  cases,  model  en¬ 
ergy  is  overpredicted  in  this  frequency  range,  limiting 
the  validity  of  the  comparison  of  the  directional 
spreading, 

•  Near  the  peak  frequency,  in  both  the  idealized  ease 
and  the  long  hindeast,  the  average  model  directional 
spreading  is  quite  dose  to  that  of  the  ground  truth. 

•  At  frequencies  above  the  spectral  peak,  the  model  in 
the  idealized  simulation  is  too  broad  relative  to  the 
ground  truth,  but  in  the  long  hindeast,  the  average 
model  directional  spreading  is  quite  close  to  that  of 
the  ground  truth. 

Regarding  statistics  other  than  bias,  in  the  long  hind- 
cast  at  all  frequency  bands,  model-data  agreement  is 
not  favorable  except  in  high-energy  conditions. 
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